/*
   S.L.U.G. : Stochastically Light Up Galaxies
  
   This code currently implements the following kludges
 	1) if a star has mass less than 0.151 Msun it is set
                to have a mass of 0.151 msun 
 	3) Stars with masses less than 0.9 (stored in variable called 
                evolved_limitmass) are not evolved. this introduces an
		approximately 0.2 dex 
               error in their measurements for approximately a Gyr
 	4) Post main sequence evolution is dealt with by a simple kill line...
               (i.e. no post main sequence evolution)
  	5) the seed is timemachine .. never tested for real stochasticity
        6) Working with fluxes in proper units, so numbers are big.
	   We should probably
           define some arbitray units to work with smaller numbers.
	7) Check what happens when timestep*sfr exceeds mass of the galaxy	
       
       
   Units:
	Bolometric luminosity: ergs/s.
	Q(H0):	s^-1
	EVERYTHING ELSE: ergs/s/Hz
   Note: throughout the code qho, qh0, and ha are all used to represent the same 
	quantity: the ionizing photons per second
*/

//include the library
#include "def.h"

using namespace std;

//global variables (set by the param file)
float maxtime;      // maximum time for which the code runs
float timestep;     // length of time steps
float galaxy_sfr;   // star formation rate in msun/yr
string galaxy_sfh;   // file name for star formation history 
string model_name;  // output model name
string imf_type;    // type of IMF to use
float metallicity; // metallicity for the stellar tracks 
string libpath;      //path to the stellar library (tree below data/)
double max_cluster_mass;   //maximum mass for the cluster mass function
                           //this doubles as the igimf_int if IGIMF is set
double slope_cmf;          //slope of the cluster mass function
double imf_slope;	   //slope for IMF if SLOPE_IMF is chosen
double min_cluster_mass;   //minimum mass for cluster mass function
double clust_frac;         //fraction of stars that are clustered
double clust_disrupt;         //parameter controlling cluster disruption
double log_time;		//this switch turns logarithmic timesteps on/off
string drawtype;   //switch to sequential/random cluster drawing1
vector<string> allfluxes;  //fluxes requested
vector<string> color1;     //color1 requested    
vector<string> color2;     //color2 requested 
vector<int> color1_ind;    //indexes for color tables     
vector<int> color2_ind;  
vector <float> allmincolor;   //min colors to use
vector <float> allmaxcolor;   //string of max colors to use
vector <float> allminflux;   //string of min fluxes to use
vector <float> allmaxflux;   //string of max fluxes to use
long bincolor;     //bin of color to use
long binflux;      //bin of flux to use
string outcluster; //enable/disable the output of the full cluster file
string imfout; //enable/disable the output of the IMF distribution of each cluster
time_t seed;    //seed to be used 

//pointer to the IMF function (set by IMF initialise)
double (*imf)(gsl_rng *) = NULL; 

//initialize some vector objects for the purpose of storing our tables
//the code works by interpolating off of a 2-dimensional tables where
//one axis is time and one is mass. Thus, given a mass and age of the
//star we can determine its properties.

//additionally we use a similar table for determining the properties
//of the low mass stars whose evolution we do not follow. These are 
//stored in variables called lowmass_****
vector<double> tab_time_axis; //the time axis of the 2d tables
vector<double> tab_time_axis_indx; //indices of the time axis
vector<double> tab_time_cutoff; //cutoff ages for main sequence stars
vector<double> tab_mass_axis; //the mass axis of the 2d tables
long nmass_axis; //the maximum size of the mass axis
vector<double> low_fluxes;
vector<double> tab_low_mass_axis;//the mass axis of the low mass stars table
long which_fluxes=0; //bitmask placeholder which keeps tracks of which fluxes 
		     //this gets set in histo_initialise.cpp!!!!


double cluster_mass=0;    //the mass of the specific cluster as determined by
                                 //the cluster mass function... this is a global variable
				// in case you use igimf_imf
double total_simulation_mass=0;

double kludge_index;
long ntables=0; //how many tables are we using
long n_flux_matrix=0;//counters for the sizes of the entries for one of the axes of these tables
long n_low_matrix=0;
//specify the parameter for the cluster mass function
string data_dir;


/*******************************************
Here the main begins
*********************************************/

int main(int argc, char *argv[]){

  double** flux_matrix;
  double** low_matrix;
  //parse the parameter file (this set the global variables!!)
  parse_param(argc,argv);

  //-----------------------------START THE JOB---------------------//
  cout<<"SLUG: Done with the global stuff!"<<endl;

  //---Some Stuff for writing out the output files---
  //make the output folder if not there
  system("mkdir -p output");
  
  //define the output filename
  string fileout("./output/");      //the output file
  fileout.append(model_name.data());
  fileout.append(".dat");
  
  string filehisto("./output/histo_");        //the output 
  filehisto.append(model_name.data());
  filehisto.append(".dat");
  
  string fileguide("./output/guide_");       //the guide for parsing the output
  fileguide.append(model_name.data());	     //the code writes out an overall output file
  fileguide.append(".dat");		     //and also this file which tells the parser 
					     //what is what

  string filetot("./output/integral_");   //Write the integrated luminosity at each time step in a smaller file
                                          //The field is also added. This is simulating the entire galaxy
  filetot.append(model_name.data());
#ifdef ASCII
  filetot.append(".txt");
#else
  filetot.append(".dat");
#endif /* ASCII */
  
  string fileimf("./output/imf_");   //write out the IMF of each cluster
  fileimf.append(model_name.data());
  fileimf.append(".dat");

  string filemisc("./output/misc_");   //write out integrated numbers that are computed at the end of the loop
  filemisc.append(model_name.data());
  filemisc.append(".dat");



  //some time machine info (also for seed)
  time_t rawtime;
  struct tm * timeinfo;
  time (&rawtime);
  timeinfo=localtime(&rawtime);  
  cout<<"SLUG: Start at time: "<<asctime(timeinfo)<<endl;
  time_t t1=time(NULL);
   
  //initialize random sequence
  const gsl_rng_type * T;
  gsl_rng * rand;
  gsl_rng_env_setup();
  if(seed == (time_t) 0) gsl_rng_default_seed=t1; else gsl_rng_default_seed=seed;
  T = gsl_rng_default;
  rand = gsl_rng_alloc (T);
  cout<<"SLUG: Seed used "<<gsl_rng_default_seed<<endl;
  
  //initialize the variables here
  double galaxy_current_mass=0;  //the current galaxies mass meaning sum of
                                 // all formed clusters
  double actual_cluster_mass=0;  //the mass one obtains by actually summing
                                 // all the masses of stars. Not always the same
				 //as galaxy_current_mass since integer number
				 //of stars is not exactly the same mass
  vector<double> star_masses;    //an array of the masses of the stars
  vector<double> star_indexes;   //the index array for the star masses

  double star_index0=0;
  double star_index1=0;
  double evolve_limitmass=0.9;  //1.5         //minimum mass of the star we evolve
  double bornage_cluster=0;              //when the cluster appears in the galaxy frame
  double evaluate_time_cluster=0;        //time at which evaluate the cluster in the galaxy frame
  double kcount=0;			//variable for holding the number of stars in a cluster
  double biggest_star=0;		//variable for holding the largest mass star in cluster
  

  double tinx_interp=0.; //interpolate time index for the star cluster
  long ntime_axis;
  double buffer,deathmass; //deathmass keeps tracks of the maximum mass star
  long lbuffer=0;
  //still alive at each timestep
  int cluster_dead=0;	//bit if cluster is disrupted or not
  long ntime;    //number of times steps
  
  double first_timestep=0;//a value used for the output guide file
  double how_many_times=0;//how many timesteps evaluated for a given cluster

  double dt_t=0; //to kill cluster during the first timestep
  long imf_type_flag=0;//this is a flag that allows us to parse which IMF
		       //we are using


  //Deal with the time step sequence and throw warnings or errors
  //set number fo time steps and the time line
  if(log_time == 0){
  ntime=(long) ceil(maxtime/timestep+1);
  } else {
  ntime=(long) ceil(log10(maxtime/log_time)/timestep +1) ; 
  }
  double *galaxy_timeline=(double *) malloc(ntime*sizeof(double));
  long nsfh_hist=(long) ceil(maxtime/1e6+1);
  double *sfh_hist=(double *) malloc(nsfh_hist*sizeof(double));
//  double *sfh_hist=(double *) malloc(10*sizeof(double));

  for(int i=0;i<nsfh_hist;i++){
    sfh_hist[i]=0;
  }

  if(maxtime>1E9){
    cout<<"SLUG: Exceeded max age 1E9 for individual star evolution and cluster distruption!!!"<<endl;
    exit(1);
  }
  
  //Do a bunch of stuff for the table and histograms
  //This is to make the code flexible to compute as many colors as you like 
  check_tables();
  cout<<"SLUG: Use "<<ntables<<" flux tables"<<endl;
  cout<<"SLUG: Compute "<<color1.size()<<" colors"<<endl;
   
  //initialize some placeholders for IMF histogramgs
  long imf_holdct;
  long imf_holdct2;
  double nimf_bin=10.0;
  double imf_hold[10];
  double imf_val[10];
  imf_val[0]=      1.46806;
  imf_val[1]=      2.39465;
  imf_val[2]=      3.90609;
  imf_val[3]=      6.37151;
  imf_val[4]=      10.3930;
  imf_val[5]=      16.9528;
  imf_val[6]=      27.6530;
  imf_val[7]=      45.1069;
  imf_val[8]=      73.5771;
  imf_val[9]=      120.017;

//  imf_val[0]=0.151;
//  imf_val[1]=0.9;
//  imf_val[2]=8;
//  imf_val[3]=13;
//  imf_val[4]=16;
//  imf_val[5]=25;
//  imf_val[6]=35;
//  imf_val[7]=50;
//  imf_val[8]=80;
//  imf_val[9]=120;
  imf_holdct=0;

  
  //Make an histogram of 10 bin in log space for the initial cluster mass function
  double cmf_binsize=0.1*log10(max_cluster_mass);
  double cmf_mass[11];     //this is for the lower limit mass in that bin. Add 1 bin for clusters
                           //that equal or are slightly bigger than max cluster mass
  double   cmf_number[11];
  int    cmf_index;
  
  //fill in the CMF mass array
  for(int i=0;i<11;i++){
    cmf_mass[i]=i*cmf_binsize;
    cmf_number[i]=0;
  }

  //initialze histograms for colors
  long nhist=(long)color1.size(); //number of colors
  
  //allocate a bunch of memory for histograms
  double * xmin=(double *) malloc(sizeof(double)*nhist);
  double * xmax=(double *) malloc(sizeof(double)*nhist);
  double * xbin=(double *) malloc(sizeof(double)*nhist);
  
  double * ymin=(double *) malloc(sizeof(double)*nhist);
  double * ymax=(double *) malloc(sizeof(double)*nhist);
  double * ybin=(double *) malloc(sizeof(double)*nhist);
  
  
  long * hf_ind1=(long *) malloc(sizeof(long)*nhist);//NOTE must be indexes of the fluxes
  //that are being used not of all fluxes
  long * hf_ind2=(long *) malloc(sizeof(long)*nhist);
  
  //set the binsize for histo
  long nx=bincolor; //nx must be the same for all histograms 
  long ny=binflux; //ny must be the same for all the histograms
  //fill with user ranges for color
  for(int kk=0;kk<(int)color1.size();kk++){
    xmin[kk]=(double)allmincolor[kk];
    xmax[kk]=(double)allmaxcolor[kk];
    hf_ind1[kk]=(long)color1_ind[kk];
    hf_ind2[kk]=(long)color2_ind[kk];
    ymin[kk]=(double)allminflux[kk];
    ymax[kk]=(double)allmaxflux[kk];
  }



  //find the "physical" x and y axis array for the color-mag diagram
  for(long histo_ct=0;histo_ct<nhist;histo_ct++){
    xbin[histo_ct]=(xmax[histo_ct]-xmin[histo_ct])/((double)nx);
    ybin[histo_ct]=(ymax[histo_ct]-ymin[histo_ct])/((double)ny);
  } 
  //now define the nasty 4d histogram vector (nx, ny, ntime, hhist)
  double **** histogram=(double ****) malloc(sizeof(double)*nx);
  for (long hct=0;hct<nx;hct++)
    histogram[hct]=(double ***) malloc (sizeof(double)*ny);
  
  for (long hct=0;hct<nx;hct++){
    for (long hct2=0;hct2<ny;hct2++){
      histogram[hct][hct2]=(double **) malloc (sizeof(double)*ntime);
    }
  }
  
  for (long hct=0;hct<nx;hct++){
    for (long hct2=0;hct2<ny;hct2++){
      for (long hct3=0;hct3<ntime;hct3++){
	histogram[hct][hct2][hct3]=(double *) malloc (sizeof(double)*nhist);
      }
    }
  }
  //initialize to zero
  for (long hct=0;hct<nx;hct++){
    for (long hct2=0;hct2<ny;hct2++){
      for (long hct3=0;hct3<ntime;hct3++){
	for (long hct4=0;hct4<nhist;hct4++){
	  histogram[hct][hct2][hct3][hct4]=0;
	}
      }
    }
  }
  
  double * big_buffer=(double *) malloc(sizeof(double)*ntables); //used later for copy/paste
  


  //load or set SFH
  //first allocate some placeholders
  double *sfh_time=NULL; // (double *) malloc(sizeof(double));
  double *sfh_cumul=NULL; //(double *) malloc(sizeof(double));
  int npntsfh=0; //how many points in SFH
  setsfh(&sfh_time,&sfh_cumul,&npntsfh);


  total_simulation_mass=sfh_cumul[npntsfh-1]; //the expected total amount of
  //mass created is just the
  //last element of the cumulative SFH
  //This is a crucial number for the code.This number 
  //controls how long we create stars. Here we turn the SFH
  //into the total number of stars created which regulates
  //how long the code runs. If you set a constant SFH, this
  //been consistently computed with the code running until
  //maxtime.


  //initialize the interpolation function for time indexes
  int intsfh=(int) npntsfh;
  gsl_interp_accel *accsfh = gsl_interp_accel_alloc ();
  gsl_spline *intrpsfh = gsl_spline_alloc (gsl_interp_linear,intsfh);
  gsl_spline_init(intrpsfh,sfh_cumul,sfh_time,intsfh);
  //to evaluate eval=gsl_spline_eval(intrpsfh,evaluate_mass,accsfh);
  
  gsl_interp_accel *accsfr = gsl_interp_accel_alloc ();
  gsl_spline *intrpsfr = gsl_spline_alloc (gsl_interp_linear,intsfh);
  gsl_spline_init(intrpsfr,sfh_time,sfh_cumul,intsfh);

 
  //prompt estimated time
  cout<<"SLUG: The estimated time to complete the program is "<<
    2*2.43*(ntime/20.)*(total_simulation_mass/4E7)*ntables/4<<" minutes"<<endl;
 

  //open the out datafile
  ofstream outfile (fileout.data(),ofstream::binary);
  ofstream outguide (fileguide.data(),ofstream::binary);
  ofstream outhist (filehisto.data(),ofstream::binary);
#ifdef ASCII
  ofstream outtot;
  outtot.open(filetot.data());
#else
  ofstream outtot (filetot.data(),ofstream::binary);
#endif /* ASCII */ 
  ofstream outimf (fileimf.data(),ofstream::binary);
  ofstream outmisc (filemisc.data(),ofstream::binary);
  
  if(imfout.compare("yes") == 0){
    //write the beginnings of the IMF file
    imf_holdct=0;
    for(int kk=0;kk<(int)nimf_bin;kk++){
      imf_hold[kk]=0;
    }
    buffer=nimf_bin+1;
    outimf.write((char *)(&buffer),sizeof(buffer));
    buffer=(double)0.9;
    outimf.write((char *)(&buffer),sizeof(buffer));
    for(int kk=0;kk<(int)nimf_bin;kk++){
      buffer=imf_val[kk];
      outimf.write((char *)(&buffer),sizeof(buffer));
    }
  }
  
  
  //write the which_fluxes bitmask
  outguide.write((char *)(&which_fluxes), sizeof(which_fluxes));
  //write out the seed
  outguide.write((char *)(&t1), sizeof(t1));
  //write out maxtime
  outguide.write((char *)(&maxtime), sizeof(maxtime));
  //write out timestep
  outguide.write((char *)(&timestep), sizeof(timestep));
  //write out the sfr
  outguide.write((char *)(&galaxy_sfr), sizeof(galaxy_sfr));


  //make a first decision on the IMF
  if (imf_type.compare("KROUPA") == 0) {
    imf_type_flag=0;
  }
  if (imf_type.compare("SALPETER") == 0) {
    imf_type_flag=1;
  }
  if (imf_type.compare("CHABRIER") == 0) {
    imf_type_flag=2;
  }
  if (imf_type.compare("IGIMF") == 0) {
    imf_type_flag=3;
  }
  if (imf_type.compare("SLOPE") == 0) {
    imf_type_flag=4;
  }


  //write out the imf_type
  outguide.write((char *)(&imf_type_flag), sizeof(imf_type_flag));

  //write out the metallicity
  outguide.write((char *)(&metallicity), sizeof(metallicity));


  //write the number of timesteps
  outfile.write((char *)(&ntime),sizeof(ntime));
#ifdef ASCII
  outtot<<ntime<<"\t"<<ntables<<"\t"<<endl;
#else
  outtot.write((char *)(&ntime),sizeof(ntime));
  //write out the number of fluxes
  outtot.write((char *)(&ntables),sizeof(ntables));
#endif
  
  //if you are fine with time (i.e. not evolving past the point where stars you
  //aren't evolving die) then go ahead and make the array of times where 
  //you want to evaluate the clusters
  if(log_time == 0){
  for(long t=0;t<ntime;t++){
    galaxy_timeline[t]=t*timestep;
    //write the number of timesteps to a file
      buffer=galaxy_timeline[t];
      outfile.write((char *)(&buffer),sizeof(buffer)); 
  }
  }else{
  for(long t=0;t<ntime;t++){
    galaxy_timeline[t]=(double) 
		log_time*pow((double)10., (double)t*timestep); //log timesteps
    //write the number of timesteps to a file
      buffer=galaxy_timeline[t];
      outfile.write((char *)(&buffer),sizeof(buffer));
  } 
  }//end log_time switch
 
  //we need to allocate memory for the flux tables which are 
  //stored in a matrix called "flux_matrix"
  flux_matrix=(double**)malloc(sizeof(double)*ntables);
  for (long i_flux_matrix=0;i_flux_matrix < ntables;i_flux_matrix++){
    flux_matrix[i_flux_matrix]=(double *) malloc (sizeof(double)*n_flux_matrix);
  } 
  low_matrix=(double**)malloc(sizeof(double)*ntables);
  for (long i_low_matrix=0;i_low_matrix < ntables;i_low_matrix++){
    low_matrix[i_low_matrix]=(double *) malloc (sizeof(double)*n_low_matrix);
  }

  //open tables with tracks and luminosity
  load_table(flux_matrix,low_matrix);
 

  double ** field_flux_matrix=(double **) malloc(sizeof(double)*ntables);
  for (long imatrix=0;imatrix<ntables;imatrix++){
    field_flux_matrix[imatrix]=(double*) malloc (sizeof(double)*ntime);
  }
  //get some sizes
  ntime_axis=(long)tab_time_axis.size();
 
  
  //now define the matrix of the fluxes
  double ** cluster_flux_matrix=(double **) malloc(sizeof(double)*ntables);
  for (long imatrix=0;imatrix<ntables;imatrix++){
    cluster_flux_matrix[imatrix]=(double*) malloc (sizeof(double)*ntime);
  }
 
  //now define the matrix of the integrated fluxes
  //this is a buffer that contains the 'galaxy'
  //integrated flux for the integrated flux output file (clusters+field)
  double ** integrated_flux_buff=(double **) malloc(sizeof(double)*ntables);
  for (long imatrix=0;imatrix<ntables;imatrix++){
    integrated_flux_buff[imatrix]=(double*) malloc (sizeof(double)*ntime);
  }
  //initialize this buffer and the field
  for(long ymatrix=0;ymatrix<ntime;ymatrix++){
    for (long xmatrix=0;xmatrix<ntables;xmatrix++){	
      integrated_flux_buff[xmatrix][ymatrix]=0.;
      field_flux_matrix[xmatrix][ymatrix]=0.;
    }
  }
  
  //initialise the pointer for the IMF
  imf_initialise(imf_type);

  //initialize the interpolation function for time indexes
  int interpoint=(int) tab_time_axis_indx.size();
  gsl_interp_accel *acc = gsl_interp_accel_alloc ();
  gsl_spline *interp = gsl_spline_alloc (gsl_interp_linear,interpoint);
  
  //the interp function doesn't accept a vector object
  double *xint_time=(double*) malloc(interpoint*sizeof(double));
  double *yint_time=(double*) malloc(interpoint*sizeof(double));
  
  for(int i=0;i<interpoint;i++){
    xint_time[i]=tab_time_axis[i];
    yint_time[i]=tab_time_axis_indx[i];
  }
  
  gsl_spline_init(interp,xint_time,yint_time,interpoint);
  

  low_fluxes.resize(ntables); //keeps track of the flux from the low mass stars
 
  kludge_index=get_tab_indexes_low(tab_low_mass_axis,0.156);

  //do stuff (untill maximum age)    
  //we now begin the loop over clusters
  cout<<"SLUG: Total Galaxy Mass= "<<total_simulation_mass<<endl;
  cout<<"SLUG: I'm computing "<<model_name<<"..."<<endl;
  
  double igimf_t0=0;
  double igimf_t1=0;
  double igimf_int=max_cluster_mass/2.;
  // interval of time to average SFR over to get an
  //estimate for the maximum cluster mass for IGIMF 
  //calculations// this uses the input for stellar
  //mass if the IGIMF is used then that input is 
  //meaningless so it doubles here
  //NOTICE WHAT I AM DOING HERE AND MAKE NOTE OF IT IN THE MANUAL


  //*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
  //*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
  //*-*-*-*-*-*-*		BEGIN ACTUAL COMPUTATION       -*-*-*-*-*-*-*-*-*
  //*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
  //*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
  //the stopping criterion is when you have built up all the stellar
  //mass that you simulation is supposed to build up

#ifndef QUIET
  float ctten=(1/31.)*total_simulation_mass;
  long  ctten0=1;
#endif /*QUIET*/ 
  
  while(galaxy_current_mass < total_simulation_mass){
#ifndef QUIET
    //Progres bar
    if(galaxy_current_mass > ctten){
      //	 cout<<ctten0<<" percent completed"<<endl;
      while(galaxy_current_mass > ctten){
	ctten0+=1;
	slug_progress(ctten0);
	ctten+=(1/31.)*total_simulation_mass;
      }
    }
#endif /* QUIET */
    
    //initialize low mass counters to 0
    for(long ntab=0;ntab<ntables;ntab++){
      low_fluxes[ntab]=0;
    }
    //now we will determine the ages of the stars in one of two ways
    // either RANDOM or SEQUENTIAL
    //here we do the 'RANDOM' message
    if (drawtype.compare("RANDOM")==0){
      bornage_cluster=gsl_spline_eval(intrpsfh, gsl_rng_uniform(rand)*sfh_cumul[npntsfh-1] ,accsfh);
    }
    //now we do the "sequential" method which is deprecated if selecteed
    if (drawtype.compare("SEQUENTIAL")==0)
      bornage_cluster=gsl_spline_eval(intrpsfh, galaxy_current_mass ,accsfh);

    if (imf_type.compare("IGIMF") == 0) {//now to deal with IGIMF
      igimf_t0=(bornage_cluster-igimf_int > 0 ? bornage_cluster-igimf_int :0);
      igimf_t1=(bornage_cluster+igimf_int < maxtime ? bornage_cluster+igimf_int: maxtime);
      galaxy_sfr=(gsl_spline_eval(intrpsfr,igimf_t1, accsfr)-gsl_spline_eval(intrpsfr,igimf_t0, accsfr))/
	(igimf_t1-igimf_t0);
      //this line calculates how many stars are created during the interval between igimf_t0 and 
      //igimf_t1 and divides by the time difference to get an instaneous SFR
      max_cluster_mass=84793.*pow((double)galaxy_sfr, (double) .75);  
      // taken from Pflam-Altenberg & Kroupa 2008 Supplementary Materials
    }

    if((1.-clust_frac)*total_simulation_mass <= galaxy_current_mass){
    //draw cluster mass
    cluster_mass=cluster_draw(rand); 
    
    //check if adding the cluster mass gets you closer to final desired mass than leaving it out
    if(fabs(total_simulation_mass-(galaxy_current_mass+cluster_mass)) >
	fabs(total_simulation_mass-galaxy_current_mass)) break;
   
    //populate the cluster
    actual_cluster_mass=cluster_populate(cluster_mass,evolve_limitmass,
					 &star_masses,&star_indexes,
					 rand, &kcount, &biggest_star,low_matrix );
 
    if(actual_cluster_mass == 0) continue; //sometimes we draw zero mass in a cluster
    if(star_masses.size() == 0) continue;

    //We use the actual cluster mass to estimate the galaxy mass
    galaxy_current_mass+=actual_cluster_mass;
    //cout<<"Drawing Cluster Mass..."<<"\t"<<actual_cluster_mass<<"\t"<<galaxy_current_mass<<"\t"<<total_simulation_mass<<endl;

    //Increment the total galaxy initial cluster mass function
    //This is taking into account all the clusters ever formed in the galaxy
  
    cmf_index=(int)floor(log10(actual_cluster_mass)/cmf_binsize);
    cmf_number[cmf_index]++;
    

    if(outcluster.compare("yes")==0){
      //write properties to the file
      outfile.write((char *)(&actual_cluster_mass),sizeof(actual_cluster_mass));//write out actual mass
      buffer=(double) kcount;//kcount holds the number of stars
      outfile.write((char *)(&buffer),sizeof(buffer));//write out number of stars
      buffer=(double) biggest_star;//this is the mass of the most massive star in the cluster
      outfile.write((char *)(&buffer),sizeof(buffer));
      buffer=(double) bornage_cluster;//this is the mass of the most massive star in the cluster
      outfile.write((char *)(&buffer),sizeof(buffer));
      //note that for the field -1 gets written for each of these quantities... so
      //if you want to add a new one be sure to add it to the final field write out...
    }
    //sort the mass array in ascending order.
    sort(star_masses.begin(),star_masses.end());  

    //initialize imf_hold to 0
     for(int kk=0;kk<(int)9;kk++){    
        imf_hold[kk]=0;
      }
     imf_holdct=0;
     
//
//begin IMF cataloging... we also use the indexing above for
//stellar evolution purposes
//
if(imfout.compare("yes") == 0){
	buffer=kcount-star_masses.size();
        outimf.write((char *)(&buffer),sizeof(buffer));
    for(int thistar=0;thistar<(int)star_masses.size();thistar++){
      if(star_masses[thistar] > imf_val[imf_holdct]){
	imf_holdct2=imf_holdct;
	while(star_masses[thistar]>imf_val[imf_holdct]){
	  imf_hold[imf_holdct]=thistar;
	  imf_holdct++;
	}
       if(imf_holdct !=0) imf_holdct--;
        imf_hold[imf_holdct2]=thistar;//no minus one because start at 0
	imf_holdct++;
	}
     }
    if(imf_holdct != (long) nimf_bin){//-1
	for(int kk3=(long)imf_holdct;kk3<(long)nimf_bin;kk3++){
          imf_hold[kk3]=star_masses.size();
	}
     }
	buffer=imf_hold[0];
        outimf.write((char *)(&buffer),sizeof(buffer));

      for(int kk=1;kk<(int) nimf_bin;kk++){
	buffer=imf_hold[kk]-imf_hold[kk-1];
        outimf.write((char *)(&buffer),sizeof(buffer));
      }
}
//
// END IMF cataloguing
//

    //time evolution
    cluster_dead=0;//a flag to denote when a cluster has been dispersed
} else {//this else is to allow non-clustered SF
    actual_cluster_mass=star_draw(cluster_mass,evolve_limitmass,
                                         &star_masses,&star_indexes,
                                         rand, &kcount, &biggest_star,low_matrix );    
    galaxy_current_mass+=actual_cluster_mass;//actual_cluster_mass in this case is now is just the stellar mass
    cluster_dead=1; // we want to add to field population

}//end else case for non-clustered star formation
    //cout<<(long) floor(bornage_cluster/1e6+0.5)<<"\t"<<bornage_cluster<<endl;
    sfh_hist[(long) floor(bornage_cluster/1e6+0.5)]+=actual_cluster_mass;
 

    first_timestep=-1; //a flag to initilaize the first time when the cluster is born
    how_many_times=0;//a counter to keep track of how many times 
		     //the cluster writes out its fluxes as a cluster...written to guide file

    for(long time=0;time<ntime;time++){//loop over timesteps for the cluster
//	cout<<"timestep"<<"\t"<<time<<endl;
      //current time for the cluster
       star_index0=(double)nmass_axis-1;
       star_index1=(double)nmass_axis-1;
      evaluate_time_cluster=galaxy_timeline[time]-bornage_cluster;
      //skip the evolution of the clusters that are not born yet
      if(evaluate_time_cluster>0){//if cluster has been born
	if(first_timestep == -1){//if cluster has been born and first_timestep is not set then set it
	  first_timestep=time;
	}
	//initialise new time step with low mass stuff
	//initialize the flux matrix with low flux values... low flux stars are NOT evolved!!!
	for (long xmatrix=0;xmatrix<ntables;xmatrix++){
	  for(long ymatrix=0;ymatrix<ntime;ymatrix++){
	    cluster_flux_matrix[xmatrix][ymatrix]=low_fluxes[xmatrix];
	  }
	}
	//get time index for the tables
	if(star_masses[0] < 1000){
	tinx_interp=gsl_spline_eval(interp,evaluate_time_cluster,acc); //this can be skipped if only lowmass

	if(tinx_interp < 0) tinx_interp=1e-3; // since tables don't run to zero time we have 
    	//cout<<"Finding time index "<<tinx_interp<<endl;
	// some problems with negative
	// indices we choose to kludge this to be set to zero 


	//find the maximum "alive" star and get rid of the dead ones
	int index;
//ISO: I had to change how this works. now we have a maximum stellar mass
// at each timestep in the tables 
//	for(index=0; index<(int)tab_time_cutoff.size()-1;index++){
//	  if(tab_time_cutoff[index] > evaluate_time_cluster) break;
//	}
	for(index=0; index<(int)tab_time_axis.size()-1;index++){
	  if(tab_time_axis[index] > evaluate_time_cluster) break;
	}
	//below is how we kill stars. if they above a designated "deathmass"
	//then they aren't evaluated. the deathmass is read from a table
	//called tab_time_cutoff...this is simply a time saving implementation
 //       if(3>4) deathmass=0.0001;
//	else 
        if(evaluate_time_cluster < tab_time_cutoff[0]) deathmass=1000.;
	//if the age of the cluster is younger than the lifetime
	//of the most massive star then the deathmass is set to 1000 msun
	//effectively turning off the deathmass cut
	else deathmass=((tab_time_cutoff[index]-tab_time_cutoff[index-1])/
			(tab_time_axis[index]-tab_time_axis[index-1]))*
	  (evaluate_time_cluster-tab_time_axis[index-1])+tab_time_cutoff[index-1];
//	else deathmass=((tab_mass_axis[index]-tab_mass_axis[index-1])/
//			(tab_time_cutoff[index]-tab_time_cutoff[index-1]))*
//	  (evaluate_time_cluster-tab_time_cutoff[index-1])+tab_mass_axis[index-1];
    
//    	cout<<"Calculating Deathamss"<<"\t"<<deathmass<<endl;
//	cout<<tab_time_axis[index-1]<<"\t"<<evaluate_time_cluster<<"\t"<<tab_time_axis[index]<<endl;
//	cout<<tab_time_cutoff[index-1]<<"\t"<<tab_time_cutoff[index]<<endl;
	//now, loop over all the remaining stars and get values
	for(int sta=0;sta<(int)star_masses.size();sta++){
	    //now get the indices. Note that star_indexes are currently 
	    //not computed in cluster_populate(). These indices are used for the
	    //the interpolation of the flux tables. Now we need 2 indices
	    // that varay as a function of timestep
//	    for(int thistar=0;thistar<(int)star_masses.size();thistar++){
//	     }
//	cout<<"mass = "<< star_masses[sta]<<"\t"<<sta<<endl;  
	  //break the for loop when you reach the mass of dead stars (this is the reason why we sorted)
//cout<<star_masses[0]<<"\t"<<deathmass<<endl;
	  if(star_masses[sta]>deathmass) break;
//cout<<"this should not print low mass"<<actual_cluster_mass<<endl;
//	      star_index0=get_star_index(tab_mass_axis,star_masses[sta], 
//				ntime_axis, floor(tinx_interp), star_index0);
//	      star_index1=get_star_index(tab_mass_axis,star_masses[sta], 
//				ntime_axis, ceil(tinx_interp), star_index1);
	//cout<<star_masses[sta]<<"\t"<<time<<"\t"<<evaluate_time_cluster<<endl;
	      star_index0=get_star_index(star_masses[sta], 
				ntime_axis, floor(tinx_interp), star_index0);
	      star_index1=get_star_index(star_masses[sta], 
				ntime_axis, ceil(tinx_interp), star_index1);

//	  cout<<"calculating star indices"<<"\t"<<star_index0<<"\t"<<star_index1<<
//		"\t"<<sta<<"\t"<<star_masses.size()<<endl;
	  for(long ntab1=0;ntab1<ntables;ntab1++){
	    big_buffer[ntab1]=interpolate_iso_phot(tinx_interp,star_index0,
  				      star_index1,ntime_axis, flux_matrix,ntab1);
	  }//big_buffer is an array with number of elements equal to ntable which has all
	  //photometry for the cluster at this timestep
	  if (nhist > 0){//if there are any histograms then do this to fill it up

	    for(long hct00=0;hct00<nhist;hct00++){
	
	      hist4d(histogram, time, 
		     big_buffer[hf_ind1[hct00]]-big_buffer[ hf_ind2[hct00]], //hf_inds store the index relative to the fluxes 
		     big_buffer[ hf_ind2[hct00]],
		     &xmin[hct00],&xmax[hct00],&xbin[hct00],&ymin[hct00],
		     &ymax[hct00],&ybin[hct00], nx, ny, hct00);//increments the histogram
	    }
	  }//end histograms
	  for(long ntab1=0;ntab1<ntables;ntab1++){
	    cluster_flux_matrix[ntab1][time]+=pow(10,big_buffer[ntab1]);//add star's flux to cluster
	  }
	}
  
      
      
      //-------------------------------------------------||
      //begin cluster death prescription
      //if the value from the uniform random number generator is not above the probability
      //of death then the cluster dies
     
  
      if(cluster_dead == 0){
	//if cluster has been born and is not dead then see if it is dead
	dt_t=timestep;	//the time interval is just the timestep unless it is the first timestep
	if(how_many_times == 0){
	  dt_t=(galaxy_timeline[time]-bornage_cluster)-1e6;
	  //this is to deal with the fact that clusters aren't alive
	  //for the entirety of their first timestep. the 1e6 is because
	  // we start killing after 1 Myr
	}

  	if((gsl_rng_uniform (rand) > clust_disrupt*dt_t/evaluate_time_cluster)
	   || (evaluate_time_cluster < 1E6)){  //start killing after 1Myr
	  //if cluster is alive then write info in the file
	  if(outcluster.compare("yes")==0){
	    for(long ntab10=0;ntab10<ntables;ntab10++){
	      buffer=cluster_flux_matrix[ntab10][time];
	      outfile.write((char *)(&buffer),sizeof(buffer));
	    }
	  }
	  how_many_times++; //keeps tracks of how many timesteps for guide file
	  

	}else{
	  cluster_dead=1; //otherwise cluster is dead
	}
      }	
      if(cluster_dead == 1){//if cluster is dead add photometry to field population
	for(long ntab1=0;ntab1<ntables;ntab1++){
	  field_flux_matrix[ntab1][time]+=cluster_flux_matrix[ntab1][time];

	}
      }
      //end cluster death prescription
      //-------------------------------------------------||
      }//end loop checking for low only

      //write the integrated flux file 
      //write all the clusters no matter if they are alive or dead to account for the field.
      for (long xmatrix=0;xmatrix<ntables;xmatrix++)
	  integrated_flux_buff[xmatrix][time]+=cluster_flux_matrix[xmatrix][time];
	
      }else{
	//fill in negative values in the arrays for times not evaluated
	for(long ntab1=0;ntab1<ntables;ntab1++){
	  cluster_flux_matrix[ntab1][time]=1E-50;
	}
      }

      
    }//end for over time
        
 
    //clear stuff that is appended before new loop  
    star_masses.clear();
    star_indexes.clear();
    
    if(outcluster.compare("yes")==0){
      //write out guide file   
      buffer=cluster_mass;
      outguide.write((char *)(&buffer),sizeof(buffer));
      
      buffer=how_many_times;
      outguide.write((char *)(&buffer),sizeof(buffer));
      
      buffer=first_timestep;
      outguide.write((char *)(&buffer),sizeof(buffer));
    }
    
  }//end for over clusters
  //add the final field to the integrated flux and write out
  for(long ymatrix=0;ymatrix<ntime;ymatrix++){
#ifdef ASCII
    outtot<<galaxy_timeline[ymatrix]<<"\t";
#endif /* ASCII */    
    for (long xmatrix=0;xmatrix<ntables;xmatrix++){	
      buffer=integrated_flux_buff[xmatrix][ymatrix];
#ifdef ASCII
      outtot<<buffer<<"\t";
#else
      outtot.write((char *)(&buffer),sizeof(buffer));
#endif /*ASCII*/
    }
#ifdef ASCII
    outtot<<endl;
#endif /*ASCII*/
  }
  
  if(outcluster.compare("yes")==0){
    //now write out the field at each timestep
    buffer=(double) -1;//write to guide file
    outguide.write((char *)(&buffer),sizeof(buffer));
    
    buffer=(double) ntime;//write to guide file
    outguide.write((char *)(&buffer),sizeof(buffer));
    
    buffer=(double) 0;//write to guide file
    outguide.write((char *)(&buffer),sizeof(buffer));
    
    buffer=(double) -1; //write to the actual file
    outfile.write((char *)(&buffer),sizeof(buffer));
    buffer=(double) -1; //write to the actual file
    outfile.write((char *)(&buffer),sizeof(buffer));
    buffer=(double) -1; //write to the actual file
    outfile.write((char *)(&buffer),sizeof(buffer));
    buffer=(double) -1; //write to the actual file
    outfile.write((char *)(&buffer),sizeof(buffer));
    
    for(long ctdum=0;ctdum<ntime;ctdum++){//loop thru timesteps for field population
      for(long ntab3=0;ntab3<ntables;ntab3++){
	buffer=field_flux_matrix[ntab3][ctdum];
	outfile.write((char *)(&buffer),sizeof(buffer));
      }
    }	
  }
  
  if(nhist != 0){
    cout<<"SLUG: Writing histograms..."<<endl;
    //here we write out the histograms
    lbuffer=(long) nx; //write out nx
    outhist.write((char *)(&lbuffer),sizeof(lbuffer));
    
    lbuffer=(long) ny; //write out ny
    outhist.write((char *)(&lbuffer),sizeof(lbuffer));
    
    
    lbuffer=(long) ntime; //write out ntime
    outhist.write((char *)(&lbuffer),sizeof(lbuffer));
    

    lbuffer=(long) nhist; //write out nhist
    outhist.write((char *)(&lbuffer),sizeof(lbuffer));

    for(long hwc=0;hwc<nhist;hwc++){
      lbuffer=(long) hf_ind1[hwc];//write index of flux 1
      outhist.write((char *)(&lbuffer),sizeof(lbuffer));

      lbuffer=(long) hf_ind2[hwc];//write index of flux 2
      outhist.write((char *)(&lbuffer),sizeof(lbuffer));
   
      //write xmin,xmax, xbin
      buffer= xmin[hwc];
      outhist.write((char *)(&buffer),sizeof(buffer));

      buffer=xmax[hwc];
      outhist.write((char *)(&buffer),sizeof(buffer));

      buffer= xbin[hwc];
      outhist.write((char *)(&buffer),sizeof(buffer));

      //write ymin,ymax, ybin
      buffer= ymin[hwc];
      outhist.write((char *)(&buffer),sizeof(buffer));

      buffer= ymax[hwc];
      outhist.write((char *)(&buffer),sizeof(buffer));

      buffer= ybin[hwc];
      outhist.write((char *)(&buffer),sizeof(buffer));
      //write out the histogram values
      for(long hwx=0;hwx<nx;hwx++){
	for(long hwy=0;hwy<ny;hwy++){
	  for (long hwt=0;hwt<ntime;hwt++){
	    buffer=(double) histogram[hwx][hwy][hwt][hwc];
	    outhist.write((char *)(&buffer),sizeof(buffer));
	  }
	}
      }

    }
  } else {
    lbuffer=(long) 0;
    outhist.write((char *)(&lbuffer),sizeof(lbuffer));

  }

  //write here the total final numbers
  double bcl;
  bcl=galaxy_current_mass;
//  outmisc.write((char *)(&galaxy_current_mass),sizeof(galaxy_current_mass));
  outmisc.write((char *)(&bcl),sizeof(bcl));
  //add the cluster mass function to the misc file

  //CMF histogram
  for(int i=0;i<11;i++){
    bcl=cmf_mass[i];
    outmisc.write((char *)(&bcl),sizeof(bcl));
    bcl=cmf_number[i];
    outmisc.write((char *)(&bcl),sizeof(bcl));
  }

  //sfh_histogram
  bcl=((double) nsfh_hist);
  outmisc.write((char *)(&bcl),sizeof(bcl));
  for(int i=0;i<nsfh_hist;i++){
  bcl=sfh_hist[i];
  outmisc.write((char *)(&bcl),sizeof(bcl));
  }

  //close the files
  outfile.close();
  outtot.close();
  outimf.close();
  outguide.close();
  outhist.close();
  outmisc.close();
  // ---the misc file has the following format---
  // 1 double: the total galaxy mass created
  // 22 doubles: entries of the CMF alternating mass/number
  // 1 double: the number of entries for the sfh histogram
  // n doubles: where n is the previous entry. each of
  //		these numbers is the SFH spaced by 1 Myr

  //at the end free random generator
  gsl_rng_free (rand);
  
  //free interp function
  gsl_spline_free(interp);
  gsl_interp_accel_free (acc);
  gsl_spline_free(intrpsfh);
  gsl_interp_accel_free (accsfh);
  gsl_spline_free(intrpsfr);
  gsl_interp_accel_free (accsfr);

  tab_time_axis.clear();
  tab_time_axis_indx.clear();
  tab_time_cutoff.clear();
  tab_mass_axis.clear();
  tab_low_mass_axis.clear();

  cout<<"SLUG: Freeing Memory"<<endl;
  //free all malloc
  free(xint_time);
  free(yint_time);
  free(sfh_hist);
  free(galaxy_timeline);
  free(sfh_time);
  free(sfh_cumul);
  free(flux_matrix);
  free(low_matrix);
  free(cluster_flux_matrix);
  free(field_flux_matrix);
  free(histogram);
  free(xmax);
  free(xmin);
  free(xbin);
  free(ymax);
  free(ymin);
  free(ybin);
  free(big_buffer);
  free(integrated_flux_buff);

  //some time info
  time (&rawtime);
  timeinfo=localtime(&rawtime);  
  cout<<"SLUG: Ended at time: "<<asctime(timeinfo)<<endl;
  //time of process
  time_t t2=time(NULL);
  cout<<"SLUG: It took "<<(t2-t1)/60.<<" minutes"<<endl;
  exit(0);
}
